////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//-------------------------------- Ebrahim Foulaadvand, 13 Oct 2012 -----------------------------------   //
//                                                                                                        //
//  The routine "DiffusionFTCS" solves the heat transfer equation using forward time centered space       //
//  scheme. Dirichlet boundary condition is used.                                                         //
//                                                                                                        //
//                                                                                                        ////                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
using namespace std;

main()
{

double k=1,tau=0.01,h,L=1,TL=10,TR=30; // k=thermal conductivity, tau=time step, h=spatial grid.

//L=system length, TL=temerature at the left boundary, TR=temerature at the right boundary.

int N=60,i,n,Time=5000; // N=number of spatial grids, Time=number of timesteps, n=timestep number.


vector<double> Tcur(N+1,0),Tnew(N+1,0); // Arrays "Tcur" and "Tnew" store the current and updated temperature values at grid points.


ofstream file1 ("Tepmerature n=10.plt");   //output file for temperature profile at timestep n=10.
ofstream file2 ("Tepmerature n=100.plt");  //output file for temperature profile at timestep n=100.
ofstream file3 ("Tepmerature n=1000.plt"); //output file for temperature profile at timestep n=1000.
ofstream file4 ("Tepmerature n=5000.plt"); //output file for temperature profile at timestep n=5000.


h=L/N;

cout<<"h= "<<h<<endl; //getch();

tau=0.5*h*h/(2*k); // adjusting the timestep with spatial grid h.

cout<<"tau= "<<tau<<endl; //getch();

Tnew[0]=TL; Tnew[N]=TR; // implementing Dirichlet boundary condition.

for(i=1;i<N;i++){Tcur[i]=0;} Tcur[0]=TL; Tcur[N]=TR; //implementing the initial condition (here delta function).

Tcur[(N/2.0)]=1/h;

// ------------------ Forward time centred space algorithm ----------------------------------

for (n=1;n<=Time;n++){

for(i=1;i<N;i++){ Tnew[i]=(k*tau/(h*h))*( Tcur[i+1]+ Tcur[i-1]-2*Tcur[i] ) + Tcur[i]; }

if (n==10) for(i=0;i<=N;i++){ file1<<-L/2.0+i*h<<"  "<<Tnew[i]<<endl; } //storing in the data file.

if (n==100) for(i=0;i<=N;i++){ file2<<-L/2.0+i*h<<"  "<<Tnew[i]<<endl; } //storing in the data file.

if (n==1000) for(i=0;i<=N;i++){ file3<<-L/2.0+i*h<<"  "<<Tnew[i]<<endl; } //storing in the data file.

if (n==5000) for(i=0;i<=N;i++){ file4<<-L/2.0+i*h<<"  "<<Tnew[i]<<endl; } //storing in the data file.


for(i=1;i<N;i++){ Tcur[i]=Tnew[i]; }

}

cout<<"end." <<endl;

getch();

}